/**
 *      triangulate.cpp
 *
 *		\mainpage triangulate
 *
 *		\section sec_intention Intention
 *
 *		Cardiac MRI images are taken to generate FEM meshes. The idea is to
 * 		create meshes for different contraction states in 2 dimenions
 * 		(cross-section through the heart). This program may be apedted to 3
 * 		dimensional problems.
 *
 * 		\section sec_desc Program description
 * 		 - 	MRI images in Dicom format need to be segmented by itksnap.
 * 			This process can be done computer assisted.
 * 		 - 	Segmentation images should be saved in gipl-format
 * 		 -	These can be read by this program, which :
 *				- 	Reads segmentation images from itk snap in gipl format and
 * 					does border interpolation to higher resolution using bsplines.
 * 					From these higher resoltuion images 2d slices are generated and a
 * 					surface triangulation is made.
 *				-	Outputs data in off format that can be read by tetgen.
 *		\section sec_compile How to compile
 *		- Dependencies:
 * 			- GSL for spline computations
 * 			- my own template library libmylib.a
 *
 *		- Compile command
 * 		g++ -Wall -I~/CPP-includes -L~/CPP-includes/lib/  triangulate.cpp -lm -lmylib -o triangulate
 *
 *		\section sec_disc Disclaimer
 *      Copyright 2009 Stefan Fruhner <stefan@itp.tu-berlin.de>
 *
 *      This program is free software; you can redistribute it and/or modify
 *      it under the terms of the GNU General Public License as published by
 *      the Free Software Foundation; either version 2 of the License, or
 *      (at your option) any later version.
 *
 *      This program is distributed in the hope that it will be useful,
 *      but WITHOUT ANY WARRANTY; without even the implied warranty of
 *      MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *      GNU General Public License for more details.
 *
 *      You should have received a copy of the GNU General Public License
 *      along with this program; if not, write to the Free Software
 *      Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
 *      MA 02110-1301, USA.
 **/

//#define DEBUG
#define PROGRAM_VERSION "0.16"
#define PROGRAM_DATE    "2011"

// To use the compute spline function in myline.h we have to define USE_EINSPLINE
//#define USE_EINSPLINE

/** \section sec_changelog Changelog
 * 	- 05.03.09 initial version
 *  - 14.03.09 new function for triangulation : triangulation of all regions is now possible
 *  - 10.05.09 pseudo 3d generation for triangle files as input (creates smesh file)
 *  - 07.09.09 selection of particular slices to generate a mesh
 *  - 13.03.11 now usage of myIniFiles, added some tuneable parameters
 *  - 13.11.11 added GSL splines instead of myMath
 *  - 22.11.11 o log messages of external programs only in verbose mode
 *             o checking the number of processors before running mkregions
 *             o convert old CARP format into new CARP format (output of mkregions)
 */

#include <iostream>
#include <limits.h>

#include <sys/stat.h> // mkdir
#include <vector>
#include <set>
#include <algorithm>
#include <mylibs/gipl.h>
#include <mylibs/mymatrix.hpp>
#include <mylibs/mystring.hpp>
#include <mylibs/point.hpp> //!< a class for points
#include <mylibs/lists.h>
#include <mylibs/maps.h>

#include <mylibs/myinifiles.hpp>
#include <mylibs/myline.hpp>
#include <mylibs/mymesh.hpp>
#include <mylibs/cmdline.hpp>
#include <mylibs/myio.hpp>
#include <MCubes.hpp>

#ifdef _OPENMP /** Compute in parallelized manner*/
	#include <omp.h>
#endif

using namespace mylibs;


#define DATATYPE unsigned short // type of GIPL data, which actually supported
#define DATATYPE_NULL_ELEMENT 0

#define USE_THRESHOLD true
#define NO_THRESHOLD false

#define is ==

#ifdef DEBUG
	#define print(x) (cout << __FILE__ << ":" << __FUNCTION__ <<"() (" <<__LINE__ << ") : " << x << endl)
#else
	#define print(x) (cout << x << endl)
#endif

// all the settings are global
double resolution 					= 		250.; 		// resolution in um
double resolution_e 				= 		250.; 		// resolution in um in low resolution regions
double scaling_factor 				= 		  1.;		// factor to scale the input coordinates
double scaling_factor_time			=  		  1.;		// factor to scale the input time
double extrude_by					= resolution;		// height after extrusion
double quality 					=  		 1.4;		// quality value needed for tetgen
double ratio_of_boundary_points 	=  		 0.5;		// ratio of boundary points  #used/ #all
set<int> high_resolution_regions;						// regions that shall be meshed in higher resolution
double ratio_edge_length_e2i		=  		10.0;		// regions can have lower resolution determined by this ratio
mystring write_basename("");
mystring write_basepath("/");
bool verbose           = false;
bool call_triangle     = false;// call triangle to create triangular mesh ?
bool call_tetgen       = false;// call tetgen to create tetrahedral mesh?
bool create_meshes     = false;// output of CARP mesh files
bool add_bath          = false;// if a bath outside shall be added
bool add_bath_auto     = true;
Point add_bath_min     = Point(); // size which shall be added to the minimum
Point add_bath_max     = Point(); // ... and to the maximum of the bounding box.
bool write_gipl_files  = false;
bool write_OOFF_files  = false;
bool write_NODE_files  = false;
bool write_MESH_files  = false;// SMESH output (tetgen)
bool write_BOUND_files = false;// output of boundary splines only
bool make_one_surface  = true; // genereate only one surface within one region
bool make_reg_surface  = false;// gen region surface with a mesh subdivided in regions
bool input_is_3d       = false;/* Input file utilizes x-y-and z, if not the
								* z-component is treated as time value. Else
								* the program computes the isosurface using the
								* marching cubes algorithm */
mystring post_processing_cmd("");


// parameters for iso-surface computation
bool isosurf                 = false; 	// flag for marching cubes iso-computation
double isovalue              = 0.;		// value used for the isosurface
int iso_smoothing_method	 = 0; 		// 0: no smoothing, 1: Laplace, 2: volume conserving LaplaceHC, 3: 1 Laplace iteration, and x LaplaceHC it
int iso_smoothing_iterations = 1; 		// Smoothness is determined by the number of iterations
float iso_laplacehc_alpha 	 = 0.1;		// parameter alpha for laplace hc smoothing algorithm
float iso_laplacehc_beta 	 = 0.3;		// parameter beta for laplace hc smoothing algorithm

bool process_only_selected_slices = false;
set<uint> slices_to_process;

//! test if an element is contained in the vector
bool contains(vector<uint> &vec, uint &value){
	for (size_t i = 0; i < vec.size(); i++){
		if (vec[i] == value) return true;
	}
	return false;
}

bool contains(set<uint> &s, uint &value){
	return (s.find(value) != s.end());
}

template <class Type, class DataType>
void make_triangulation(SurfaceMesh &mesh, Type &matrix, DataType first_value, uint z, uint t, bool use_threshold, DataType *threshold);
myList<PointCurve> do_border_interpolation(SurfaceMesh &mesh, double resolution = 0.25, double ratio = 0.2);
bool gen_border(myList<PointCurve> &Spline, int z=0, int t=0, bool only_init=false);
void make_FEM_mesh(SurfaceMesh &mesh, uint k);
void make_FEM_regions(SurfaceMesh &mesh, GIPL &img, uint k);
void make_pseudo3d(SurfaceMesh &triangulation, double z_value);
void delete_node_files(const char *base, bool verbose=false);
void convert_node_to_CARP(const char *basename);
void convert_oldCARP_to_newCARP(const char *basename);
void compute_iso_surface(SurfaceMesh &mesh, mymatrix<DATATYPE> &matrix, const double isovalue);
void extrude_mesh(myStringList &non_opt_args);

int main(int argc, char** argv){

	//settings needed to run the program
	DATATYPE threshold[2] = {0,0};//initialze with some values
	vector<DATATYPE> vthreshold;

	GIPL img;
	myIniFiles ini("triangulate.ini",argc, argv, false);

	ini.dependency_checking(true);
	ini.set_info("\n\
 >>> triangulate " + mystring(PROGRAM_VERSION) + " <<< \n\
     (compiled on " + toString(__DATE__) + " " + toString(__TIME__) + ")\n\
 Generates a surface triangulation of 2D gipl images.\n\
 - Triangulates the surface of an object coded in the gipl file by a\n\
   certain color range  (= threshold)\n\
 - 3d and 4d gipl files are interpreted as slices in z direction\n\
   with time values. A mesh is created for each z and t !\n\
 - saves the surface off the object as off-format\n\
 - creates pseudo 3d data by copying each slice to z-resolution\n\
 - calls tetgen to create a volume tetrahedralization\n\
\n					(c) Stefan Fruhner "+mystring(PROGRAM_DATE)+"\n\
 usage:\n\
	../triangulation <infile>.gipl output [other params (see below)]\n\
	\n\
	<infile>.gipl is the file to be read and output is the \n\
	basenane of the mesh to be written, e.g. folder/mesh. The folder\n\
	will be created.\n\
	Attention: The data type of the GIPL is must be 'unsigned short'.");

	ini.register_flag ("verbose"     , "v", "Detailed console messages");
	ini.register_param("rescale"     , "c", "<float> - rescale all coordinates of the input file by the factor, e.g. input is given in mm, but output is needed in um. (default = 1.) \n\nAttention: the resolution must be defined in the units after the scaling factor is applied !");
	ini.register_param("threshold"   , "t", "<int>,<int> : treshold color values which numerically enclose the object to be triangulated") ;
	ini.register_param("edge_ratio"  , "e", "<float> - ratio of resolution in high resolution regions to low resolution regions (default = 10.) ");
	ini.register_param("resolution"  , "r", "<float> - (high) resolution value of the output mesh");
	ini.register_param("high_res_reg", "m", " <int>  - identifier for high_resolution regions (default = 1) ");
	ini.register_param("rescale-time", "i", "<float> - rescale the time defined in the input file by a factor (default = 1.)");
	ini.register_param("ratio_bnd"   , "u", "<float> - ratio of boundary points to use for spline computation (default=0.5)");

	ini.register_param("slice"       , "s", "<int>,...  : process only slices with index given in the list");
	ini.register_flag ("extrude"     , "X", "<filename> : generate (P)seudo 3d from 2d-mesh files by extruding existing 2D surface (for example from triangle .node / .ele files). Overrides all params except -x");
	ini.register_param("extrude_by"  , "x", "<float> : height in z-direction after extrusion (after rescaling).");
	ini.register_flag ("regions"     , "R", "Create a triangulation with regions (-b or -C mustn't be set)");
	ini.register_flag ("bath"        , "B", "Emerge the extracted object in a surrounding bath. To the boundary box a certain distance is added. By default this is 0., but your can specify own values in bath_add_min and bath_add_max.");
	ini.register_param("bath_add_min", "a" , "Size to add to the minium  extension of the object (default = auto). Implies -B.");
	ini.register_param("bath_add_max", "A" , "Size to add to the maxmimum extension of the object (default = auto). Implies -B.");
	ini.register_flag ("boundaries"  , "b", "Save spline interpolation of border only. Saves .node files containing all boundary curves. Each curve gets an own label. Overrides most of the other flags (e.g. -gCIMNROST)");
	ini.register_flag ("CARP"        , "C", "As -b but additionally creates meshes from boundary curves -- writes CARP-type files.");
	ini.register_flag ("NODE"        , "N", "By default meshes are created in Carp-format. You can additionally keep the node-files created by triangle and tetgen.");
	ini.register_flag ("GIPL"        , "g", "Save *.gipl file with the defined resolution");
	ini.register_flag ("OFF"         , "O", "Save *.off file (written out by tetgen). Contains all elements. (geomview format)");
	ini.register_flag ("SMESH"       , "S", "Save *.smesh-file. Contains surface only (with region labels). Can be used for fast visualization. (It is also used as tetgen input) ");
	ini.register_flag ("triangle"    , "I", "call triangle to create triangular FEM mesh from boundaries.");
	ini.register_flag ("tetgen"      , "T", "call tetgen to create tetrahedral FEM mesh from extruded surface.") ;
	ini.register_param("quality"     , "q", "<float> - quality of the mesh (edge to surface ratio for tetgen)", 1.4 , 0.2, 10., "tetgen") ;
	ini.register_param("post_processing", "p", "Shell script containing commands to post process the mesh files Occurrences of %%meshbase%% and %%logfile%% will be replaced apropriately.");

	ini.register_flag("2D", "2", "Input data is only two dimensional with z as time evolution");

	ini.register_flag("isosurf", 0, "Computes the isosurface using the marching cubes algorithm. Input date must be three-dimensional.");
	ini.register_param("isovalue", 0, "Value used to compute iso-surface.", true, "isosurf");
	ini.register_param("smoothing_method", 0, "Method used for surface smoothing \n [0 = 'no smoothing', 1 = 'Laplace', 2 = 'volume conserving  LaplaceHC', 3 = '1xLaplace+n*LaplaceHC'].", 0, 0, 3, "isosurf");
	ini.register_param("smoothing_iterations", 0, "Number of smoothing iterations. The higher the value the smoother the meshes.",true, "isosurf");
	ini.register_param("laplacehc_alpha", 0, "Parameter alpha for LaplaceHC smoothing algorithm.", 0.1, 0., 1., "smoothing_method>1");
	ini.register_param("laplacehc_beta", 0, "Parameter beta for LaplaceHC smoothing algorithm.", 0.3, 0., 1., "smoothing_method>1");

	ini.check(2);

	ini.read("threshold", vthreshold);
	if (vthreshold.size() == 2 ){
		//! \note Error checking is done later, we have to distinguish between -P or not
//		cmdline::exit("Error: Threshold level must have 2 values (-t <lower bound>,<upper bound>)");
		threshold[0] = vthreshold[0];
		threshold[1] = vthreshold[1];
	}

	resolution               = ini.read("resolution"  ,resolution				);
	ratio_edge_length_e2i	 = ini.read("edge_ratio"  ,ratio_edge_length_e2i	);
	scaling_factor			 = ini.read("rescale"     ,scaling_factor			);
	scaling_factor_time		 = ini.read("rescale-time",scaling_factor_time		);
	quality                  = ini.read("quality"	  ,quality					);
	ratio_of_boundary_points = ini.read("ratio_bnd"	  ,ratio_of_boundary_points	);

	// read the regions that should be meshes in high res
	high_resolution_regions.clear();
	ini.read("high_res_reg", high_resolution_regions);

	if ((ini.exists("edge_ratio") and high_resolution_regions.size() == 0)){
		cmdline::exit("You need to define either a high\
 resolution region (high_res_reg) or skip 'edge_ratio'");
	} else if ( (not (ini.exists("edge_ratio"))
					and high_resolution_regions.size() == 0)){
		cmdline::warning("Attention: No high resolution region defined !");
	}

	verbose					 = ini.exists("verbose"		);
	call_tetgen				 = ini.exists("tetgen"		);
	call_triangle			 = ini.exists("triangle"	);
	create_meshes			 = ini.exists("CARP"		);
	write_gipl_files		 = ini.exists("GIPL"		);
	write_OOFF_files		 = ini.exists("OFF"			);
	write_MESH_files		 = ini.exists("SMESH"		);
	write_NODE_files		 = ini.exists("NODE"		);
	write_BOUND_files		 = ini.exists("boundaries"	);
	add_bath				 = ini.exists("bath"		);

	if ((ini.exists("bath_add_min") or ini.exists("bath_add_max"))
			and not ini.exists("bath")){
		add_bath = true;
	}

	vector<double> abv;
	ini.read("bath_add_min", abv);
	add_bath_min.set(abv);
	ini.read("bath_add_max", abv);
	add_bath_max.set(abv);

	add_bath_auto = (not ini.exists("bath_add_min") and not ini.exists("bath_add_max"));

	post_processing_cmd 	 = ini.read("post_processing", mystring(""));
	if (not post_processing_cmd.empty() and not mylibs::IO::file_exists(post_processing_cmd.c_str()))
		cmdline::exit("post_processing_cmd file not found !");

	input_is_3d				 = not ini.exists("2D");
	// marching cubes parameters
	isosurf			 		 = ini.exists("isosurf");
	isovalue				 = ini.read("isovalue", isovalue);
	iso_smoothing_iterations = ini.read("smoothing_iterations", 1);
	iso_smoothing_method	 = ini.read("smoothing_method", 0);
	iso_laplacehc_alpha		 = ini.read("laplacehc_alpha", iso_laplacehc_alpha);
	iso_laplacehc_beta		 = ini.read("laplacehc_beta", iso_laplacehc_beta);

	// set high resolution resolution
	resolution_e = resolution * ratio_edge_length_e2i;

	if (ini.exists("regions")){
		make_reg_surface =  true;
		make_one_surface =  false;
	}

	slices_to_process.clear();
	ini.read("slice", slices_to_process);
	if (slices_to_process.size() > 0 ) process_only_selected_slices = true;

	extrude_by = ini.read("extrude_by", resolution); // attention valid for all modes

// --- parsing parameters is ready ---

	mystring non_opt_args = ini.param_non_opt_args();

	cout << "remains: "<< non_opt_args << endl;
	myStringList opts = non_opt_args.split();

	if (ini.exists("extrude")){
		extrude_mesh(opts);
		exit(EXIT_SUCCESS);
	}

	mystring file("");
	bool file_defined = false;
	for (size_t i = 0; i < opts.size(); i++){
		mystring &opt = opts[i];

		if (not file_defined){		//! first  arg = input file
			file = opt;
			file_defined = true;
		} else{ 						//! second arg = output file
			write_basename = opt;
			write_basepath = write_basename.file_path() + '/';
			mkdir(write_basepath.c_str(), 0740);
		}

	}

	if (not file_defined) cmdline::exit("\aError:	No input file given ... ");
	if (write_basename =="") write_basename = mystring("output");

	if (file.file_ext() == "node"){
		cmdline::msg(" - Loading boundaries from node file.");
		cmdline::msg("Currently not implemented");

		if (verbose){
			cout << " - Reading           file  : "	<< file 				 << endl;
			cout << " - Basename for output is  : " << write_basename		 << endl;
			cout << " - resolution              : "	<< resolution 			 << endl;
			cout << " - scaling_factor          : "	<< scaling_factor		 << endl;
			cout << " - scaling_factor_time     : "	<< scaling_factor_time	 << endl;
			cout << " - ratio edge len. ext/int : "	<< ratio_edge_length_e2i << endl;
			cout << " - spline bnd ratio        : "	<< ratio_of_boundary_points		<< endl;
			cout << " - Mesh quality            : " << quality 				 << endl;
			cout << " - extrude_by              : " << extrude_by 			 << endl;
		}
		myList<PointCurve> Boundaries;
		PointCurve c(file);
		PointCurve::split_multiple_curves(Boundaries, c);
		gen_border(Boundaries);
		exit(EXIT_SUCCESS);
	} else if (file.file_ext() == "gipl"){
		img = GIPL(file.c_str()); // read in the header of the gipl file

		// rescale coordinates
		img.origin(	img.origin(0)*scaling_factor,
					img.origin(1)*scaling_factor,
					img.origin(2)*scaling_factor,
					img.origin(3)*scaling_factor_time);

		img.pixdim(	img.pixdim(0)*scaling_factor,
					img.pixdim(1)*scaling_factor,
					img.pixdim(2)*scaling_factor,
					img.pixdim(3)*scaling_factor_time );
	} else 	cmdline::exit("\aError:	No known file type given ... ");


	if ((threshold[0] == 0) and (threshold[1] == 0))
			cmdline::exit("\aError:	Threshold levels must be defined !");

	if (verbose){
		cout << " - Reading gipl - file     : "	<< file 				 << endl;
		cout << " - Basename for output is  : " << write_basename		 << endl;
		cout << " - Using threshold values  : ["<< threshold[0]			 << ","<<  threshold[1] << "]" << endl;
		cout << " - resolution              : "	<< resolution 			 << endl;
		cout << " - scaling_factor          : "	<< scaling_factor		 << endl;
		cout << " - scaling_factor_time     : "	<< scaling_factor_time	 << endl;
		cout << " - ratio edge len. ext/int : "	<< ratio_edge_length_e2i << endl;
		cout << " - spline bnd ratio        : "	<< ratio_of_boundary_points		<< endl;
		cout << " - Mesh quality            : " << quality 				 << endl;
		cout << " - extrude_by              : " << extrude_by 			 << endl;
		img.print_header();
	}

	if ( img.image_type() != GIPL_U_SHORT )
		cmdline::exit("Other data type then unsigned short currently not supported,");

	// read all image data
	mymatrix<DATATYPE> matrix(img.dims(0),img.dims(1),img.dims(2),img.dims(3));
	cout << " - Reading data";
	matrix.read_data(file, img.header_size()); // read all image data
	cout << "   ...done." << endl;

	matrix.origin(	img.origin(0),img.origin(1), img.origin(2),img.origin(3) );
	matrix.pixdim(	img.pixdim(0),img.pixdim(1), img.pixdim(2),img.pixdim(3) );


	SurfaceMesh mesh;

	if (isosurf){ 	// if computation of iso-surfaces
		compute_iso_surface(mesh, matrix, isovalue);

		switch (iso_smoothing_method){
			case 0:	break;
			case 1:	mesh.SmoothSurfaceLaplacian(iso_smoothing_iterations); 	break;
			case 2:	mesh.SmoothSurfaceLaplacianHC(iso_smoothing_iterations, iso_laplacehc_alpha, iso_laplacehc_beta);break;
			case 3:
				mesh.SmoothSurfaceLaplacian(1);
				mesh.SmoothSurfaceLaplacianHC(	iso_smoothing_iterations,
												iso_laplacehc_beta,
												iso_laplacehc_beta);
				break;
		}

		mesh.save_mesh("test_iso");
		mesh.save_smesh("test_iso.smesh");
//					mesh.save_poly("test_iso.poly");
//					mesh.save
		exit(EXIT_SUCCESS);
	} // else continue here


	for (uint t = 0; t < matrix.dims(3); t++){			// for all times (4th dimension)
		for (uint k = 0; k < matrix.dims(2); k++){		// for all z values (3rd dimension)

			//! checks if the user wants this slice to be processed
			if (process_only_selected_slices ){
				switch (matrix.dim()){
				case 3:
					//! \attention gen_border needs the first mesh as
					//! initialization to set the border dimensions correctly
					if (k==0 and (not contains(slices_to_process,k) ))
						if (create_meshes or write_BOUND_files){
							mesh.clear(); // clean up in every slice
							print(" - Computing surface triangulation");
							make_triangulation(mesh, matrix, matrix[0], k, t, USE_THRESHOLD, threshold);
							mesh.smooth_borders2D();
							myList<PointCurve> Spline = do_border_interpolation(mesh, resolution, ratio_of_boundary_points);
							gen_border(Spline, k, t, true);
							continue;
						}

					if (not contains(slices_to_process,k) ) continue; // skip this one
					break;
				case 4:
					if (not contains(slices_to_process,t) ) continue; // skip this one
					break;
//							default:
				}
			}

			print ("\n######### [ z = " << k << ", t = " << t <<"] ######################");
			mesh.clear(); // clean up in every slice
			print(" - Computing surface triangulation");
			make_triangulation(mesh, matrix, matrix[0], k, t, USE_THRESHOLD, threshold);

			mesh.smooth_borders2D();
			Point min;
			Point max;
			mesh.bounding_box(min, max);
			min.info("min");
			max.info("max");

			if (create_meshes or write_BOUND_files){
				myList<PointCurve> Spline = do_border_interpolation(mesh, resolution, ratio_of_boundary_points);
				gen_border(Spline, k, t);
				continue;
			}
			// do the interpolation some smoothing save the
			// surface triangulation and finally call tetgen
			if (make_one_surface)	make_FEM_mesh(mesh, k);
			else
			if (make_reg_surface)	make_FEM_regions(mesh, img, k);
		}
	}

	print("ready.");
	return 0;
}

myList<PointCurve> do_border_interpolation(SurfaceMesh &mesh, double resolution, double ratio){
	print(" - Detecting border of the surface");

#ifndef __GSL_SPLINE_H__
	print(" - ALARM : You have to compile mylibs using the GSL spline library !");
	exit(EXIT_FAILURE);
#endif


	myList<PointCurve> border = mesh.find_border();	//! compute the border sorted by arc-length

	print(" - Doing a spline interpolation ( ratio = "<< ratio << " , resolution = "<< resolution << " )" );
	myList<PointCurve> Spline;
	for (size_t i = 0; i < border.length(); i++){ //! for every curve do the interpolation
		//PointCurve Interpolation = border[i].compute_spline(ratio, resolution, false); // the last bool means: saving files spline.dat, curve.dat
		PointCurve Interpolation = border[i].compute_spline(ratio, resolution,0); // the last bool means: saving files spline.dat, curve.dat
//		PointCurve Interpolation = border[i].compute_spline(1., resolution, 0, true); // the last bool means: saving files spline.dat, curve.dat
		Spline.append(Interpolation);			  //! ... and save this is one list of points
	}
	return Spline;
}

void make_FEM_mesh(SurfaceMesh &mesh, uint k){
try {
	myList<PointCurve> splines = do_border_interpolation(mesh, resolution/2., ratio_of_boundary_points);

	// merge the result into one spline
	PointCurve Spline;
	for (size_t i = 0; i < splines.size(); i++) Spline.merge(splines[i]);

	Spline.sort(); //! \attention Here we destroy the order by arc-length!;

	// Computing interpolated image
	print(" - Computing interpolated image");
	// find the minimum and maximum dimensions of the curve
	Point min = Spline.min();
	Point max = Spline.max();

	cout << "   - Top left corner at     \t:\t "; min.prnt();
	cout << "   - Bottom right corner at \t:\t "; max.prnt();

	// compute high resolution image
	double width  = max.x - min.x;
	double height = max.y - min.y;

	int xitems = int(width /resolution) + 1;
	int yitems = int(height/resolution) + 1;

	print("   - width : " << xitems << " , height : " << yitems);

	mymatrix<unsigned char> high_res(xitems,yitems);
	// set the new origin
	high_res.origin(min);
	high_res.pixdim(resolution, resolution);

	//! project the curve to a regular grid
	Point p = Spline[0];			// set the starting point for iteration
	for (size_t i = 0; i < Spline.length(); i++){
			try{
				uint index = high_res.CoordsToIndex(p);
				high_res[index] = 250;	// set some arbitrary color value
				p = Spline.next();		// iterate to the next point of the curve
			}
			catch (string e) { print(i << " " << e);} // if the indices don't match print a warning !
	}
	//compute center of mass
	print(" - Computing center of mass");
	Point com = Spline.center_of_mass();

	/** \attention Since we sorted the points we will find
	in the first item the point which is the `lexicograpically smallest`**/
	Point diff =  com - Spline[0];
	diff.normalize(); 			// we are only interested in the direction
	diff *= resolution * 1.41; 		// length should be one resolution unit times sqrt(2) (maximum distance to the next pixel);

	Point start = Spline[0] + diff;
	uint starter = high_res.CoordsToIndex(start);

	print(" - Filling the surface with color");
	high_res.flood_fill(250, starter);

//! \todo Check if the user wants to extrude the mesh

	print(" - Creating pseudo 3D images");
	high_res.append(high_res);   			// append the matrix to itself
	high_res.dims(xitems, yitems, 2); 		// we now have two slices
	high_res.pixdim(resolution, resolution, resolution , 0.);


	GIPL high_res_gipl;
	high_res_gipl.image_type(GIPL_U_CHAR); //! \attention Don't forget to set the right datatype
	high_res_gipl.dims(xitems, yitems, 2, 1); 	// we now have two slices
	high_res_gipl.pixdim(resolution, resolution, resolution , 0.);
	high_res_gipl.origin(min);

	high_res.origin(min);
	high_res.pixdim(resolution, resolution, resolution , 0.);

	//save in a file
	if (write_gipl_files) {
		char fn[100];
		sprintf(fn, "%s%02d.gipl",write_basename.c_str(), k);
		print(" - Saving to file " << fn <<".");
		FILE *out = high_res_gipl.save_header(string(fn));
		high_res.save_to_file(out);
	}

	print(" - Computing triangulation of high resolution segmentation");
	SurfaceMesh high_res_mesh;

	unsigned char thresh[2] = {100, 255}; // triangulate item between 100 and 255
	make_triangulation(high_res_mesh,high_res, high_res[0], 0, 0, USE_THRESHOLD, thresh);

	print(" - Smoothing edges");
	high_res_mesh.smooth_borders2D();

	print(" - Creating pseudo 3d suface mesh");
	high_res_mesh.pseudo3d(extrude_by, resolution);
	for (size_t i = 0; i < high_res_mesh.region_def.size(); i++){ // update region definitions
		/**! \todo Should not be applied anymore since pseudo3d centers the object during extrusion
					Check if it works ... **/
//		high_res_mesh.region_def[i].pos.z -= extrude_by / 2.;  // reduce z coordinate by half of the extrusion height
		// 5 tets fit into one cube
		if (high_resolution_regions.find(i) != high_resolution_regions.end())
			high_res_mesh.region_def[i].constraint = (resolution   * resolution   * resolution  ) / 5.;
		else
			high_res_mesh.region_def[i].constraint = (resolution_e * resolution_e * resolution_e) / 5.;
	}

	char fn[256];
	if (write_OOFF_files){
		sprintf(fn, "%s_%02d.off",write_basename.c_str(), k);
		high_res_mesh.save_off(fn);	 //! save triangulation
	}

	if (call_tetgen or write_MESH_files){
		sprintf(fn, "%s_%02d.smesh",write_basename.c_str(), k);
		high_res_mesh.save_smesh(fn);	 //! save triangulation
	}

	if (not call_tetgen) {
		// \todo save Carp meshes here
		return;
	}

	char cmd[512];
	print(" - Calling tetgen");
	if (write_OOFF_files)	sprintf(cmd,"tetgen %s -q%f -C -O", fn, quality);
	else					sprintf(cmd,"tetgen %s -q%f -C   ", fn, quality);

	cmdline::run_command(cmd); // call tetgen

}
catch (string(e)) {print(e); }
}

template <class Type>
void paint_curves(myList<PointCurve> &Curves, mymatrix<Type> &matrix){
	for (size_t j = 0; j < Curves.size(); j++){
		PointCurve &spline = Curves[j];
		Point p = spline[0];			// set the starting point for iteration

		//! project the curve to a regular grid
		for (size_t i = 0; i < spline.length(); i++){
			try{
				uint index = matrix.CoordsToIndex(p);
				matrix[index] =j+1;		// set some arbitrary color value
				p = spline.next();		// iterate to the next point of the curve
			}
			catch (string e) { print(i << " " << e);} // if the indices don't match print a warning !
		}
	}
}


void make_FEM_regions(SurfaceMesh &mesh, GIPL &img, uint k){

	try {
	myList<PointCurve> Spline = do_border_interpolation(mesh, resolution/2., ratio_of_boundary_points);

  // Computing interpolated image
	print(" - Computing interpolated image");
	// find the minimum and maximum dimensions of the curve
	Point min = img.coords_top_left();
	Point max = img.coords_bottom_right();
	cout << "   - Top left corner at     \t:\t "; min.prnt();
	cout << "   - Bottom right corner at \t:\t "; max.prnt();

	// compute high resolution image
	double width  = max.x - min.x;
	double height = max.y - min.y;

	int xitems = int(width /resolution) + 1;
	int yitems = int(height/resolution) + 1;

	print("   - width : " << xitems << " , height : " << yitems);

	mymatrix<unsigned char> high_res(xitems,yitems);
	// set the new origin
	high_res.origin(min);
	high_res.pixdim(resolution, resolution);

	SurfaceMesh high_res_mesh;
	GIPL high_res_gipl;
	high_res_gipl.image_type(GIPL_U_CHAR); //! \attention Don't forget to set the right datatype
	high_res_gipl.dims(xitems, yitems, 1, 1); 	// we now have two slices
	high_res_gipl.pixdim(resolution, resolution, 0., 0.);
	high_res_gipl.origin(min);

	//! Draw curve on regular grid
	paint_curves(Spline, high_res);

	//! compute for every curve a triangulation and append theses triangles
	//! to the mesh
	for (size_t j = 0; j < Spline.size(); j++){
		PointCurve spline = Spline[j];

		spline.sort();
			//compute center of mass
		print(" - Computing center of mass");
		Point com = spline.center_of_mass();
		spline.sort(); 						//! sort the curve lexicographically
		/** \attention Since we sorted the points we will find
		in the first item the point which is the `lexicograpically smallest`**/
		Point diff =  com - spline[0];
		diff.normalize(); 			// we are only interested in the direction
		diff *= resolution * 1.41; 		// length should be one resolution unit times sqrt(2) (maximum distance to the next pixel);

		size_t start = high_res.CoordsToIndex(spline[0] + diff);

		print(" - Filling the surface with color");
		high_res.flood_fill(j+1, start);
	}

	unsigned char thresh[2] = {0,0};
	make_triangulation(high_res_mesh, high_res, high_res[0], 0, 0, NO_THRESHOLD, thresh);

	high_res_mesh.scan_for_regions();

//! \todo check if the user wants tu extrude the mesh, otherwise save different formats here ad exit.

	print(" - Creating pseudo 3D images");
	high_res.append(high_res);   			// append the matrix to itself
	high_res.dims(xitems, yitems, 2); 		// we now have two slices
	high_res.pixdim(resolution, resolution, resolution , 0.);

	//save to a file
	if (write_gipl_files) {
		high_res_gipl.dims(xitems, yitems, 2); 		// we now have two slices
		high_res_gipl.pixdim(resolution, resolution, resolution , 0.);

		char fn[100];
		sprintf(fn, "%s%02d.gipl",write_basename.c_str(), k);
		print(" - Saving to file " << fn <<".");
		FILE *out = high_res_gipl.save_header(string(fn));
		high_res.save_to_file(out);
	}

	print(" - Smoothing edges");
	high_res_mesh.smooth_borders2D();

	print(" - Creating pseudo 3d suface mesh");
	high_res_mesh.pseudo3d(extrude_by, resolution);
	for (size_t i = 0; i < high_res_mesh.region_def.size(); i++){
		/**! \todo Should not be applied anymore since pseudo3d centers the object during extrusion
					Check if it works ... **/
//		high_res_mesh.region_def[i].pos.z -= extrude_by / 2.;
		// 5 tets fit into one cube
		if (high_resolution_regions.find(i) != high_resolution_regions.end())
			high_res_mesh.region_def[i].constraint = (resolution   * resolution   * resolution  ) / 5.;
		else
			high_res_mesh.region_def[i].constraint = (resolution_e * resolution_e * resolution_e) / 5.;
	}

	char fn[256];
	if ( write_OOFF_files){
		sprintf(fn, "%s_%02d.off",write_basename.c_str(), k);
		high_res_mesh.save_off(fn);	 //! save triangulation with colored regions
	}

	if (call_tetgen or write_MESH_files){
		sprintf(fn, "%s_%02d.smesh",write_basename.c_str(), k);
		high_res_mesh.save_smesh(fn);
	}

	if (not call_tetgen) return;

	char cmd[512];
	if (write_OOFF_files)	sprintf(cmd,"tetgen %s -q%f -C -O", fn, quality);
	else					sprintf(cmd,"tetgen %s -q%f -C   ", fn, quality);

	cmdline::run_command(cmd); // call tetgen
	}
	catch (string(e)) {print(e); }
}

/** name: gen_border
 * Determines boundaries and computes spline interpolations, saves them
 * as *.poly file, calls triangle, tetgen and does optional post
 * processing in order to create meshes.
 *
 * @param mesh : mesh of the segmentation's surface
 * @param z : z-coordinate of the slice
 * @param t : time value of the slice
 * @param only_init : static_values will be initialized only
 * @return bool (false on error)
 **/
bool gen_border(myList<PointCurve> &Spline, int z, int t, bool only_init){
//	myList<PointCurve> Spline = do_border_interpolation(mesh, resolution, ratio_of_boundary_points);

// create a poly file for triangle

	//! append outer boundary as curve
	PointCurve outer_bnd;
	static  Point top_left  = Point(0.,0.,0.);  		//! static means here: only set in the first run, so that d is eqaul for all meshes
	static  Point bot_right = Point(0.,0.,0.);
	static double zval = 0.;
	static bool first_call = true;

	Spline[0].min().info("Spline min");
	Spline[0].max().info("Spline max");

	if (first_call){ // initializiations on first function call
		first_call = false;
		zval = Spline[0].max().z; 							//! z-value of mesh is taken from the zeroth z-slice
		if (add_bath_auto){
			Point d   = Spline[0].max() - Spline[0].min();
			top_left  = Spline[0].min() - d/2.; 			//! edges of the mesh	min - d/2
			bot_right = Spline[0].max() + d/2.; 			//!							max + d/2
			if (not input_is_3d){
				top_left.z  = zval;
				bot_right.z = zval;
			}
		}else {
			top_left  = Spline[0].min() - add_bath_min; 	//! edges of the mesh	min
			bot_right = Spline[0].max() + add_bath_max; 	//!							max
		}

	}

	if (only_init) return true;

	myList<int> nodecount; // count nodes per region
	//	save boundary curves to file
	{
		PointCurve curves;
		// determine the number of nodes in each region
		for (size_t region = 0; region < Spline.size(); region++){
			PointCurve &spline = Spline[region];
			spline.attribute((int)region);
			int cnt =  ( spline.circular() ) ? (int)spline.size()-1 : (int) spline.size() ;
			curves.extend(spline);
			nodecount.append( cnt);
		}

		char fn_node[512];
		sprintf(fn_node, "%s_%02d_%02d_curves.node",write_basename.c_str(),  z, t);

		// set fixed z-value (original z-data has the meaning of time)
		if (not input_is_3d)
			for (size_t i = 0; i < curves.size(); i++) curves[i].z = zval;

		curves.save_node(fn_node, 3);
	}
	//! \todo Split function into upper part (computation of boundary) and lower part (further processing).
	if (write_BOUND_files) return true; // skip the rest of the function
	// add a bath
	if ( add_bath == true){
		Point top_right(top_left);	top_right.x = bot_right.x;
		Point bot_left(bot_right);	bot_left.x  = top_left.x;

		//append corners
		outer_bnd.append(top_left);
		outer_bnd.append(top_right);
		outer_bnd.append(bot_right);
		outer_bnd.append(bot_left);
		outer_bnd.append(top_left);

		if (verbose){
			if (not input_is_3d)
				cout << " - Setting z value  = "    << zval << endl;
			cout << " - setting outer borders " << endl;
			top_left.info("top left    ");
			top_right.info("top right   ");
			bot_right.info("bottom right");
			bot_left.info("bottom left ");
		}

		if (not call_tetgen or ( add_bath_min.z == add_bath_max.z and add_bath_min.z == 0. )){
			Spline.push_front(outer_bnd);
			nodecount.push_front(4);
		}
	}

	int nodes = 0;
	for (int i = 0; i < (int)nodecount.size(); i++)	nodes += nodecount[i];

	char fn[512];
	sprintf(fn, "%s_%02d_%02d.poly",write_basename.c_str(),  z, t);
	print(" - Writing *.poly file : " << fn);

	FILE *ou = 0;
	ou = fopen(fn,"w");
	if (!ou)return false;

//! *  First line: <# of vertices> <dimension (must be 2)>  <# of attributes> <# of boundary markers (0 or 1)>
	fprintf(ou, "%d 2 0 1\n", nodes);

//! * Following lines: <vertex #> <x> <y> [attributes] [boundary marker]
	int offset = 0;
	for (size_t region = 0; region < Spline.size(); region++){
		Point &pt = Spline[region][0]; // select the first point
		for (int i = 0; i < nodecount[region]; i++){
				fprintf(ou, "%d %f %f %d\n", offset + i, pt.x, pt.y, (int) region+1 );
				pt = Spline[region].next();
			}
		offset += nodecount[region];
	}

	int segments = 0;
	for (int i = 0; i < (int)nodecount.size(); i++)	{
		segments += nodecount[i];
		if (not Spline[i].circular()) segments--;
	}
	fprintf(ou, "# line segments with boundary markers\n");
//!    * One line: <# of segments> <# of boundary markers (0 or 1)>
	fprintf(ou, "%d 1\n", segments);

//!    * Following lines: <segment #> <endpoint> <endpoint> [boundary marker]
	offset = 0;
	for (size_t region = 0; region < Spline.size(); region++){
		for (int pt = 0; pt < nodecount[region]-1; pt++)
				fprintf(ou, "%d %d %d %d\n", offset + pt , offset + pt, offset + pt + 1, (int)region+1);

		if ( Spline[region].circular() )
			fprintf(ou, "%d %d %d %d\n", offset + nodecount[region] -1, offset + nodecount[region]-1, offset, (int)region+1);
		else
			offset--; // if the curve is not circular

		offset += nodecount[region];
	}

//! * One line: <# of holes>
//! * Following lines: <hole #> <x> <y>
	fprintf(ou, "# no holes\n");
	fprintf(ou, "0\n");

//! * Optional line: <# of regional attributes and/or area constraints>
//! * Optional following lines: <region #> <x> <y> <attribute> <maximum area>
	fprintf(ou, "# regions \n");
	fprintf(ou, "%d\n", (int)Spline.size());

	print(" - Setting regions");
	fprintf(ou, "# <region #> <x> <y> <attribute> <maximum area> \n");
	for (int i = 0; i < (int)Spline.size() ; i++){
		Spline[i].sort();

		// select point inside the regions
		Point com = Spline[i].center_of_mass();
		/** \attention Since we sorted the points the first point is
		 * `lexicograpically smallest`**/
		Point diff =  com - Spline[i][0];
		diff.normalize(); 				// we are only interested in the direction
		diff *= resolution * 1.41; 		// length should be one resolution unit times sqrt(2) (maximum distance to the next pixel);

		Point start = Spline[i][0] + diff;
		if (not input_is_3d) start.z = zval;
		double area_constraint = 0.5 * resolution_e * resolution_e;		 // area constraint, for low res regions

		// high resolution region
		if (high_resolution_regions.find(i) != high_resolution_regions.end())
			area_constraint = 0.5 * resolution * resolution; // 0.5 * (edge length)^2

		fprintf(ou, "%d %f %f %d %f \n", i, start.x, start.y, i	, area_constraint );

		printf("   Region %d area constraint = %.6f at ", i, area_constraint); start.prnt();
	}
	fprintf(ou, "\n");
	fclose(ou);

	if (not call_triangle) return true;

	char cmd[512];
	if (verbose)
		sprintf(cmd, "triangle -pq30.VAa %s_%02d_%02d.poly > %s_%02d_%02d_triangle.log",write_basename.c_str(), z, t,write_basename.c_str(), z, t);
	else
		sprintf(cmd, "triangle -pq30.VAa %s_%02d_%02d.poly > /dev/null",write_basename.c_str(), z, t);
//	sprintf(cmd, "triangle -pq25.VAa %s_%02d_%02d.poly > %s_%02d_%02d_triangle.log",write_basename.c_str(), z, t,write_basename.c_str(), z, t);
	cmdline::run_command(cmd);


	//! create a pseudo 3d mesh
	// read in again from triangle files
	SurfaceMesh triangulation;
	char base[1024];
	sprintf(base,"%s_%02d_%02d.1",write_basename.c_str(),z,t);

	triangulation.read_mesh(base);

	if (not call_tetgen){ // if only 2d meshes shall be created
		sprintf(base,"%s_%02d_%02d",write_basename.c_str(),z,t);
		triangulation.save_mesh(base);
	}

//	make_pseudo3d(triangulation, top_left.z);
	make_pseudo3d(triangulation, zval);


	if (not write_NODE_files) { // remove all triangle files
		sprintf(base,"%s_%02d_%02d.1",write_basename.c_str(),z,t);
		delete_node_files(base, verbose);
		sprintf(base,"%s_%02d_%02d.1.poly",write_basename.c_str(),z,t);
		mylibs::IO::rm(base);
	}

	sprintf(base,"%s_%02d_%02d.poly",write_basename.c_str(),z,t);
	mylibs::IO::rm(base);

	//! \todo incorporate bath in z-direction, if demanded
	if ( add_bath == true and not ( add_bath_min.z == add_bath_max.z and add_bath_min.z == 0 )  ){
		// relabel all regions, since we want to append the bath as region 0
		for (size_t i = triangulation.region_def.size(); i > 0; i--){
			triangulation.exchange_face_attribute(i-1, i);
			triangulation.region_def[i-1].attribute = i;
		}

		Point p0(bot_right); p0.x = top_left.x;		p0.boundary(1);
		Point p1(bot_right); 						p1.boundary(1);
		Point p2(bot_right); p2.y = top_left.y;		p2.boundary(1);
		Point p3(top_left ); p3.z = bot_right.z;	p3.boundary(1);

		Point p4(p0); p4.z = top_left.z;			p4.boundary(1);
		Point p5(p1); p5.z = top_left.z;			p5.boundary(1);
		Point p6(p2); p6.z = top_left.z;			p6.boundary(1);
		Point p7(top_left);							p7.boundary(1);

		//append corners
		int i0 = triangulation.points(); triangulation.append_point(p0);
		int i1 = triangulation.points(); triangulation.append_point(p1);
		int i2 = triangulation.points(); triangulation.append_point(p2);
		int i3 = triangulation.points(); triangulation.append_point(p3);
		int i4 = triangulation.points(); triangulation.append_point(p4);
		int i5 = triangulation.points(); triangulation.append_point(p5);
		int i6 = triangulation.points(); triangulation.append_point(p6);
		int i7 = triangulation.points(); triangulation.append_point(p7);

		triangulation.append_face(Face(i0,i1,i3,0));
		triangulation.append_face(Face(i1,i2,i3,0));
		triangulation.append_face(Face(i1,i5,i2,0));
		triangulation.append_face(Face(i5,i6,i2,0));

		triangulation.append_face(Face(i4,i5,i6,0));
		triangulation.append_face(Face(i4,i6,i7,0));
		triangulation.append_face(Face(i4,i7,i3,0));
		triangulation.append_face(Face(i0,i4,i3,0));

		triangulation.append_face(Face(i0,i1,i4,0));
		triangulation.append_face(Face(i1,i5,i4,0));
		triangulation.append_face(Face(i3,i2,i7,0));
		triangulation.append_face(Face(i2,i6,i7,0));

//		triangulation.append_region(SurfaceMesh::Region(top_left, 4, 3200000000.000000));
		triangulation.region_def.push_front(SurfaceMesh::Region(top_left, 0, 3200000000.000000));

	}

	if (write_OOFF_files){
		sprintf(base,"%s_pseudo3d_%02d_%02d.off",write_basename.c_str(),z,t);
		triangulation.save_off(base);
	}

	sprintf(base,"%s_3d_%02d_%02d.smesh",write_basename.c_str(),z,t);
	triangulation.save_smesh(base);

	if (not call_tetgen) return true;

	if (write_OOFF_files){ // to preserve surface, add Y : No splitting of input boundaries (facets and segments)
		if (verbose)
			sprintf(cmd, "tetgen -OpqA %s > %s_%02d_%02d_tetgen.log", base,base, z, t);
		else
			sprintf(cmd, "tetgen -OpqA %s > /dev/null", base);
	}
	else{
		if (verbose)
			sprintf(cmd, "tetgen -pqA %s > %s_%02d_%02d_tetgen.log", base,base, z, t);
		else
			sprintf(cmd, "tetgen -pqA %s > /dev/null", base);
	}
	cmdline::run_command(cmd);

	// delete the smesh file
	if (not write_MESH_files) mylibs::IO::rm(base);

	// post processing step
	// input mesh ==> tetgen style
	if (not post_processing_cmd.empty()){
		if (not mylibs::IO::file_exists(post_processing_cmd.c_str()))
			cmdline::exit("post_processing_cmd file not found !");

		char cm[1024];
		ifstream pp(post_processing_cmd.c_str(), ios::in);
		mystring cmd("");
		while (pp.good()){
			pp.getline(cm, 1024);
			cmd += mystring(cm)+"\n";
		}
		pp.close();

		// run
		sprintf(base,"%s_3d_%02d_%02d",write_basename.c_str(),z,t);
		cmd.replace("%%meshbase%%", base);

		mystring meshdir(base);
		cmd.replace("%%meshdir%%", meshdir.file_path().c_str());

		if (cmd.contains("%%logfile%%")){
			char logfile[1024];
			sprintf(logfile, "%s_%02d_%02d_mkregions.log", base , z, t);
			cmd.replace("%%logfile%%", logfile);
		}
		cmdline::run_command(cmd.c_str());
//		//!
//		sprintf(base,"%s_3d_%02d_%02d_t",write_basename.c_str(),z,t);
//		convert_oldCARP_to_newCARP(base);

//		sprintf(base,"%s_3d_%02d_%02d_t_i",write_basename.c_str(),z,t);
//		convert_oldCARP_to_newCARP(base);
	}
//	if (call_mkregions){
//		//! \todo use mkregions_cmd_file
//		// run mkregions
//		sprintf(base,"%s_3d_%02d_%02d",write_basename.c_str(),z,t);
//		if (verbose)
//			sprintf(cmd, "mpirun -np %d `which mkregions` %s > %s_%02d_%02d_mkregions.log  < /dev/null ", (int) num_procs, base, base , z, t);
//		else
//			sprintf(cmd, "mpirun -np %d `which mkregions` %s > /dev/null  < /dev/null ",(int) num_procs, base);

//		cmdline::run_command(cmd);
//		sprintf(base,"%s_3d_%02d_%02d_t",write_basename.c_str(),z,t);
//		convert_oldCARP_to_newCARP(base);

//		sprintf(base,"%s_3d_%02d_%02d_t_i",write_basename.c_str(),z,t);
//		convert_oldCARP_to_newCARP(base);
//	}
	else{ // then we have node files here
		sprintf(base,"%s_3d_%02d_%02d.1",write_basename.c_str(),z,t);
		convert_node_to_CARP(base);
	}

	if (not write_NODE_files) {
		sprintf(base,"%s_3d_%02d_%02d.1",write_basename.c_str(),z,t);
		if (not write_NODE_files) delete_node_files(base);
	}

	return true;
}

//! check if a value is in between threshold values
#define in_range(v) ( ((v >= threshold[0]) && (v <= threshold[1])) ? true:false )

/** name: make_triangulation
 * Performs a marching square algorithm to compute a triangulation of a 2D
 * region.
 * @param mesh          : the resulting triangular mesh is written to this object.
 * @param matrix        : Object containing regular data points with values with up to 4D data [x,y,z,t]
 * @param first_value   : dummy parameter that must be of the same type as the data
 * @param z             : Defines the z-value of the slice to be extracted
 * @param t             : Defines the t-value of the slice to be extracted
 * @param use_threshold : true if some threshold values are present
 * @param threshold     : values defining the region to be extracted
 * @return
 **/
template <class Type, class DataType>
void make_triangulation(SurfaceMesh &mesh, Type &matrix, DataType first_value, uint z, uint t, bool use_threshold, DataType *threshold){
	mesh.clear(); // erase mesh contents

	first_value = first_value;

	for (uint j = 0; j < matrix.dims(1)-1; j++) 		//! y - loop
	for (uint i = 0; i < matrix.dims(0)-1; i++){		//! x - loop
		size_t p[4];
		p[0] = matrix.index(i,j,z,t); // index in der Matrix
		p[1] = p[0] + matrix.d(1);
		p[2] = p[1] + matrix.d(0);
		p[3] = p[0] + matrix.d(0);

		// get all values at all points
		DataType v[4];
		for (int k = 0; k < 4; k++)	v[k] = matrix[p[k]];

		/** neighbouring pixels
		 *  x 3 x
		 *  0 X 2
		 *  x 1 x
		 */
		size_t n[4];
		n[0] = (i > 0				)? p[0] - matrix.d(0) : 0;
		n[1] = (j < matrix.dims(1)-2)? p[0] + matrix.d(1) : 0;
		n[2] = (i < matrix.dims(0)-2)? p[0] + matrix.d(0) : 0;
		n[3] = (j > 0				)? p[0] - matrix.d(1) : 0;
		//! ... and their values
		DataType nv[4];
		for (int k = 0; k < 4; k++)	nv[k] = matrix[n[k]];


		bool has_neighbour[4];
		has_neighbour[0] = (i > 0				 )? true : false;
		has_neighbour[1] = (j < matrix.dims(1)-2 )? true : false;
		has_neighbour[2] = (i < matrix.dims(0)-2 )? true : false;
		has_neighbour[3] = (j > 0				 )? true : false;

		bool n_equal[2]; // connections between the edges
		n_equal[0] = (has_neighbour[3] and has_neighbour[0] and nv[3] == nv[0] ) ? true: false;
		n_equal[1] = (has_neighbour[0] and has_neighbour[1] and nv[0] == nv[1] ) ? true: false;

		if (( use_threshold == true)
		and ( not in_range(v[0] ) )
		and	( not in_range(nv[0]) )
		and ( not in_range(nv[1]) ) )
			continue;

		Point pt0 = matrix.coords(i  ,j  ,z,t);
		Point pt1 = matrix.coords(i  ,j+1,z,t);
		Point pt2 = matrix.coords(i+1,j+1,z,t);
		Point pt3 = matrix.coords(i+1,j  ,z,t);

		// append points to the mesh, edge: index in the point list
		int edge0 = mesh.append_point(pt0, p[0]);
		int edge1 = mesh.append_point(pt1, p[1]);
		int edge2 = mesh.append_point(pt2, p[2]);
		int edge3 = mesh.append_point(pt3, p[3]);


		/** append triangles
		  *   a       b        c      d
		  * X . x	X --x	X . x	X --x
		  * | \ .	. \ |	. / |	| / x 	X is the point p0
		  * x-- x	x . x	x --x	x . x	Points are numbered counter clockwise
		  */
		if ((n_equal[0]) && (v[0] != nv[0]) ) {
			if (use_threshold) {
				if (in_range(nv[0])) mesh.append_face(Face(edge0,edge1,edge3, nv[0])); //d
				if (in_range(v[0] )) mesh.append_face(Face(edge1,edge2,edge3, v[0] )); //c
			}
			else {
				mesh.append_face(Face(edge0,edge1,edge3, nv[0])); //d
				mesh.append_face(Face(edge1,edge2,edge3, v[0] )); //c
			}
			continue;
		}

		if ((n_equal[1]) && (v[0] != nv[1])) {
			if (use_threshold) {
				if (in_range(nv[1])) mesh.append_face(Face(edge0,edge1,edge2, nv[1])); //a
				if (in_range(v[0] )) mesh.append_face(Face(edge0,edge2,edge3, v[0]) ); //b
			}
			else {
				mesh.append_face(Face(edge0,edge1,edge2, nv[1])); //a
				mesh.append_face(Face(edge0,edge2,edge3, v[0]) ); //b
			}
			continue;
		}

		// if nothing matches than use the own value
		if (	(! use_threshold) or
				((use_threshold) && (in_range(v[0] )))
			){
			if ((i+j) % 2 == 0){
				// if nothing matches than use the own value
				mesh.append_face(Face(edge0,edge1,edge3, v[0]));
				mesh.append_face(Face(edge1,edge2,edge3, v[0]));
			}
			else{
				mesh.append_face(Face(edge0,edge1,edge2, v[0]));
				mesh.append_face(Face(edge0,edge2,edge3, v[0]));
			}
		}
	}

	mesh.determine_dimension(); // should be 2D
}

void make_pseudo3d(SurfaceMesh &triangulation, double z_value){

	if (not input_is_3d){
		// set the z coordinate for all points
		for (size_t i = 0; i < triangulation.points(); i++){
			triangulation.p[i].z = z_value; //top_left.z;
		}

		// set the z-coordinate for all regions
		for (size_t i = 0; i < triangulation.region_def.size(); i++){
			triangulation.region_def[i].pos.z = z_value;// top_left.z;
		}
	}

	// finally make it pseudo 3d
	triangulation.pseudo3d(extrude_by, resolution);

	for (size_t i = 0; i < triangulation.region_def.size(); i++){
		/**! \todo Should not be applied anymore since pseudo3d centers the object during extrusion
					Check if it works ... **/
//		triangulation.region_def[i].pos.z -= extrude_by / 2.;
		// 5 tets fit into one cube
		if (high_resolution_regions.find(i) != high_resolution_regions.end())
			triangulation.region_def[i].constraint = (resolution   * resolution   * resolution  ) / 5.;
		else
			triangulation.region_def[i].constraint = (resolution_e * resolution_e * resolution_e) / 5.;
	}
}

void delete_node_files(const char *base, bool verbose){
	char fn[1024];

	strncpy(fn, base, 1024);
	strncat(fn, ".ele", 1024);
	if (verbose) cmdline::msg("   Deleting file " +string(fn) + ".");
	mylibs::IO::rm(fn);

	strncpy(fn, base, 1024);
	strncat(fn, ".face", 1024);
	if (verbose) cmdline::msg("   Deleting file " +string(fn) + ".");
	mylibs::IO::rm(fn);

	strncpy(fn, base, 1024);
	strncat(fn, ".node", 1024);
	if (verbose) cmdline::msg("   Deleting file " +string(fn) + ".");
	mylibs::IO::rm(fn);

}

void convert_node_to_CARP(const char *basename){
	char base[1024];
	sprintf(base,"%s.node",basename);
	SurfaceMesh final_mesh(base);
	final_mesh.save_mesh(basename);
}

void convert_oldCARP_to_newCARP(const char *basename){
	char base[1024];
	sprintf(base,"%s.pts",basename);
	SurfaceMesh final_mesh( base, (verbose == true) );
	final_mesh.save_mesh(basename);

	sprintf(base,"%s.tetras",basename);
	mylibs::IO::rm(base);

	// copy old *.ttlon file to new *.lon file
	sprintf(base,"%s.lon",basename);
	fstream out(base, fstream::out);

	sprintf(base,"%s.ttlon",basename);
	fstream in(base, fstream::in);

	out << "1" << endl;
	out << in.rdbuf();

	in.close();
	out.close();

	// remove original ttlon file
	mylibs::IO::rm(base);
}

void compute_iso_surface(SurfaceMesh &mesh, mymatrix<DATATYPE> &matrix, const double isovalue){
	cout << " - Computing iso-surface" << endl;
	if (matrix.dim() < 3) throw myexception("Image dimension < 3");

	bool endcaps = true;
	size_t add = 0;
	if (endcaps) add = 2;
	cout << matrix.dims(0) * matrix.dims(1) * matrix.dims(2)  << " items" << endl;
	double *data  = new double[matrix.dims(0) * matrix.dims(1) * (matrix.dims(2)+add)];
	if (not data) throw myexception("Insufficent memory");

	for (uint z = 0; z < matrix.dims(2)+add; z++){
		for (uint y = 0; y < matrix.dims(1); y++){
			for (uint x = 0; x < matrix.dims(0); x++){
				uint idx2 = matrix.index(x,y,z);
				uint idx = matrix.index(x,y,z);
				if (endcaps){
					idx = matrix.index(x,y,z+1);
					if ((z == 0) or (z == matrix.dims(2)+1)) {
						data[idx2] = 0.;
						continue;
					}
				}
				data[idx2] = (double) matrix[idx];
			}
		}
	}

	mesh.clear();
	int att = 1;

	MCubes::Grid mcubes(matrix.pixdim(), matrix.dims(0), matrix.dims(1), matrix.dims(2));
	mcubes.RunAlgorithm(&data[0], isovalue);
	delete[] data; data = 0;
	list<MCubes::Triangle>::iterator it = mcubes.TriangleList.begin();
	for (it = mcubes.TriangleList.begin(); it != mcubes.TriangleList.end(); it++) {
		MCubes::Triangle &t = *it;
		size_t a = mesh.append_point_uniquely(Point( t.v1.x, t.v1.y, t.v1.z,0.,att ));
		size_t b = mesh.append_point_uniquely(Point( t.v2.x, t.v2.y, t.v2.z,0.,att ));
		size_t c = mesh.append_point_uniquely(Point( t.v3.x, t.v3.y, t.v3.z,0.,att ));
		Face f = Face(a, b, c); f.bnd = 1; f.attribute = att;
		mesh.append_face(f);
	}


	for (size_t i = 0; i < mesh.points(); i++)	mesh.p[i].boundary(1);

}


void extrude_mesh(myStringList &opts){
	print(" - Pseudo 3d computation switched ON ... ");
	string file("");
	for (size_t i = 0; i < opts.size(); i++){
		mystring &opt = opts[i];

		if (file.empty()) file = opt;
		else
		if ( write_basename.empty() ){
			write_basename = opt;					// output file
			write_basepath = write_basename.file_path() + '/';
			mkdir(write_basepath.c_str(), 0740);
		}
	}

	mystring basename = mystring(file);
	basename = basename.file_base();
	print(" - Reading files with base " << basename);
	print(" - Extrusion height        " << extrude_by);
	SurfaceMesh triangulation;
	if (write_basename.empty()) write_basename = basename+"_extruded";
	triangulation.read_mesh(basename.c_str());
	make_pseudo3d(triangulation, 0.);
	triangulation.save_smesh( (write_basename + ".smesh").c_str());
}
